1 Antes de empezar, conozcámonos

Hola, soy Aitor Gonzalez. Yo aprendo tidymodels para que te sea más fácil de aprender.

Soy data scientist o lo que es lo mismo, estadístico que hace más cosas aparte de una regresión logística. Actualmente trabajo en la unidad de facultad de medicina y biología de la Universidad Autónoma de Barcelona.

Este taller contiene mucho texto. Y suele requerir de un par de veces mirárselo para pillarle todo el tranquillo. No pasa nada si no te enteras a la primera. “Dentro de poco” estará subido a Youtube en mi canal Reestimando.

2 ¿Por que tidymodels?

2.1 Un poquito de la historia

La historia comienza con este señor de la foto de la ziquierda. Él es Max Khun. Parece un señor muy majo o muy travieso. Yo personalmente creo más lo segundo viendo su creación.

Comienza por aquí por que este señor fue el desarrollador de la librería caret. caret fue una librería que comenzó en 2008 y que asentó las bases del modelado estadístico en R. Desde su creación, el paquete ha sido citado en más de 7000 publicaciones académicas, teniendo un enorme éxito en la comunidad.

El paquete fue un éxito debido a que antes de él, no había una forma de unificar sintaxis entre los principales paquetes de modelado. Cada modelo estaba en un paquete y todos diferentes entre si, un caos para organizar un flujo de trabajo. caret permitía utilizar con una sintaxis común los modelos de varias librerías. Además también traía varias funciones para hacer automáticamente cosas que son necesarias en todo proceso de modelado, como una función que haga una partición de datos en un conjunto de entrenamiento y uno de prueba (training y test).

ELIPSIS en 2012 aparece la el entorno tidyverse de la mano de Hadley Whickam, la otra súper estrella de R en los últimos años. Tidyverse proporciona un montón de herramientas para poder lidiar con facilidad con la manipulación de datos, convirtiéndose en un package súper popular no solo ne R, sino en toda la ciencia datos. El principal punto es la facilidad de sintaxis para todas sus funciones, haciendo que los usuarios tuvieran una lectura de código muy fácil y fluida, enfocándose 100% en la necesidad y no en la programación.

Encantado con un framework en el que toda la sintaxis está ordenada, en 2019, Max se une al tidyverso, y comienza a crear tidymodels cooperando con el equipo de Whickam. El objetivo era convertir Careta algo más modularizado y que permitiera centrar en el modelado en si mismo, más que infinitas líneas de sintaxis para poder crear modelos de machine learning.

A fecha de este taller, noviembre de 2023, hay otra desarrolladora principal de tidymodels además de Max, Julia Silge.

Esta mujer, no solo es un as de la programación, sino que además se dedica a expandir la filosofía tidy a través de sus redes. Su contenido vale oro y diría que es imperdible si quieres no solo aprender más de tidymodels, sino de las posibilidades del aprendizaje automático en general.

- https://www.youtube.com/@JuliaSilge

- https://juliasilge.com/blog/

Entre algunos de sus trabajos más llamativos están uno que me encanta que es Topic modeling for #TidyTuesday Taylor Swift lyrics es decir, “modelizando los temas que tratan las canciones de Taylor Swift” https://juliasilge.com/blog/taylor-swift/

2.2 Un breve ejemplo

Imaginemos que queremos hacer una regresión lineal. Cogeremos unos datos cualquiera, lo de mtcars.supongamos que queremos hacer una regresión logítica, la sintaxis para eso sería:

Como se pude ver en el ejemplo, la regresión fuera de tidymodels la hacemos con la librería glm. En tidymodels esto se especifica con la función set_engine.

Ahora, veamos que pasaría si queremos hacerla con las liberría h2o

Se va pillando la idea, ¿no?…

tidymodels permite que con cambios mínimos de código puedas cambiar completamente tu flujo de trabajo. Si queremos cambiar un “motor” de modelo, solo tenemos que especificarlo en un solo cambio, y no tener que reescribir todo el código de nuevo. Esto es clave ya que permite jugar mucho a la hora de probar muchos tipos de modelos y ampliar las miras sin aumentar los problemas de código.

Evidentemente, esto no es ni la superficie de tidymodels, aquí solo estamos viendo la isla a lo lejos.

2.3 ¿Que aporta de tidymodels?

No más caret: caret está en mantenimiento solo, en el futuro se dará soporte a tidymodels. tidymodels es la apuesta al largo plazo en el entorno de R.
Consistencia: Misma sintaxis para todos los procesos.
Replicabilidad: Es muy sencillo replicar resultados.
Comunicación: Los outputs de tidymodels siguen la lógica de tidyverse.
Exportabilidad: Junto a otras librerías, como plumber y vetiver, tidymodels permite fácilmente la exportabilidad de sus modelos a otros entornos en cloud, para poner los modelos a producción.
Posibilidades: a fecha de este informe, tidymodels da soporte a 155 modelos, incluyendo modelos para datos censurados y modelos de series temporales.

3 Presentando los datos

Somos unos expertos en aprendizaje automático y nos piden ayuda de un hospital 🏥

Al parecer, nuestra variable, Var_respuesta, es una enfermedad complicada de diagnosticar. Los médicos no son capaces de determinar si un paciente tendrá o no esa enfermedad al cabo de un año. Entonces han hecho un estudio caso-control en el que se miran unas cuantas variables; y al cabo de un año se da un diagnostico definitivo de var_respuesta. En este caso pues, nos interesa una modelización de tipo clasificatoria binaria, es decir de 2 niveles, o tiene la enfermedad o no la tiene. 👍/👎

Además de de nuestra variable Var_respuesta, tenemos otras 25 variables de diversos tipos, de cognición, bioquímicas, genéticas, etc. El dataframe es una simulación de uno real, e incluye valores perdidos y otras características a tener en cuenta en la fase de pre-procesamiento.

Insisto, el propósito es modelizar un clasifiación binaria

Antes de empezar podemos ver un poco su estructura.

#---------------
# Librerias | Carga de datos ----
#---------------
if (!require('pacman') ) {install.packages('pacman')}
pacman::p_load(
  readxl,      # Lectura de los datos
  tidyverse ,  # Acceso al entorno de procesamiento "tidyverse"
  tidymodels,  # Acceso al entorno de modelado "tidymodels"
  agua,        # Modelaje de entorno h2o 
  themis,      # Soporte de la librería "recipes" para tratar desbalanceo
  skimr,       # Descriptiva rápida de datos  
  naniar,      # Resumen y Visualización de datos faltantes
  knitr,       # esta en html para este formato
)

#--- Carga de datos
datos <- read_xlsx('Data/datos.xlsx') %>% 
  filter(!is.na(Var_respuesta) ) %>% 
  mutate(Var_respuesta = as_factor(Var_respuesta))

#--- Primeras 10 filas de los datos.
datos %>% 
  head(.,10) %>% 
  knitr::kable()

Var_respuesta	Edad	Cognicion_01	Cognicion_02	Cognicion_03	Cognicion_04	Cognicion_05	Bioquimica_01	Bioquimica_02	Bioquimica_03	Bioquimica_04	Bioquimica_05	Bioquimica_06	Bioquimica_07	Bioquimica_08	Bioquimica_09	Bioquimica_10	Bioquimica_11	Escala_01	Escala_02	Escala_03	Genetica_01	Genetica_02	Genetica_03	Genetica_04	Genetica_05	Genetica_06	Genetica_07
0	14	101.82	NA	NA	NA	NA	NA	NA	5.65	13.65	125.71	NA	NA	84.29	1.48	30.35	35.68	0.00000	28.397013	7.363183	NA	0.9665678	NA	NA	1.2921746	0.2417377	NA
0	17	80.41	85.53	55.68	124.52	159.70	62.37	NA	5.65	17.83	154.53	NA	NA	99.42	1.34	33.99	33.01	138.73119	16.769985	17.764831	0.4550571	NA	NA	NA	NA	NA	0.0044324
0	14	91.00	79.76	94.28	140.68	100.07	86.33	26.24	3.18	13.09	135.88	81.53	5.75	109.45	0.95	30.44	34.21	NA	21.382220	7.061147	NA	NA	NA	NA	NA	NA	NA
0	21	86.99	91.33	93.65	115.00	123.01	66.84	10.23	3.48	11.97	132.15	90.66	2.58	68.60	1.04	30.10	34.00	48.83394	7.785552	20.442332	NA	NA	NA	NA	NA	NA	NA
1	33	78.83	104.19	42.68	112.97	179.16	59.51	NA	4.05	13.99	118.33	114.46	3.03	62.87	1.14	22.08	33.50	33.37826	28.165315	28.817018	NA	1.0690024	0.0339304	0.7820054	0.3230982	NA	0.3998781
1	27	93.54	81.16	89.32	115.83	65.07	93.17	32.35	3.50	17.95	132.61	67.61	4.10	76.30	1.11	36.50	32.24	33.45837	7.376139	22.080631	1.3081069	0.5763763	NA	0.5438819	NA	0.4705209	NA
0	16	114.25	NA	NA	NA	NA	NA	8.14	5.68	12.25	149.89	NA	NA	91.22	1.28	35.86	40.11	NA	18.489371	10.365972	1.5581016	NA	NA	1.1597700	NA	0.6242926	0.0579575
0	13	74.80	NA	NA	NA	NA	NA	93.32	4.77	8.88	164.33	NA	4.33	84.49	1.07	23.16	27.76	NA	18.749210	9.068513	NA	NA	NA	NA	NA	NA	NA
0	14	66.78	NA	NA	NA	NA	NA	54.33	2.56	16.07	161.64	NA	NA	99.47	0.92	30.14	34.27	133.14928	18.116769	31.680105	1.6247538	NA	0.9775707	NA	0.9699578	NA	1.1983120
1	17	84.13	100.23	69.62	137.84	91.24	78.76	28.80	4.22	14.38	151.67	NA	5.23	91.19	1.75	18.24	37.76	113.35402	20.327740	21.101902	NA	NA	1.1801013	NA	NA	1.0568568	NA

#--- Valoración general
skimr::skim(datos)

Data summary
Name	datos
Number of rows	285
Number of columns	28
_______________________
Column type frequency:
factor	1
numeric	27
________________________
Group variables	None

Variable type: factor

skim_variable	n_missing	complete_rate	ordered	n_unique	top_counts
Var_respuesta	0	1	FALSE	2	0: 227, 1: 58

Variable type: numeric

skim_variable	n_missing	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
Edad	9	0.97	22.33	6.47	8.00	17.00	21.00	27.00	39.00	▂▇▆▅▂
Cognicion_01	0	1.00	86.51	29.37	0.00	78.50	90.56	103.13	143.84	▁▁▅▇▂
Cognicion_02	82	0.71	90.22	12.98	55.23	80.74	90.43	98.98	129.15	▁▆▇▃▁
Cognicion_03	57	0.80	77.62	17.80	35.22	65.91	76.77	88.23	130.09	▂▇▇▃▁
Cognicion_04	64	0.78	134.85	51.16	15.31	97.17	137.91	172.32	245.38	▂▅▇▆▂
Cognicion_05	70	0.75	126.23	44.38	53.89	94.13	120.96	155.32	249.33	▆▇▆▂▁
Bioquimica_01	54	0.81	75.94	14.24	42.80	66.92	74.56	84.50	118.52	▂▇▇▃▁
Bioquimica_02	43	0.85	36.07	38.35	0.12	11.43	24.47	42.21	250.61	▇▂▁▁▁
Bioquimica_03	13	0.95	4.05	1.64	1.27	2.90	3.67	4.74	12.68	▇▇▂▁▁
Bioquimica_04	12	0.96	14.30	1.96	8.88	12.95	14.34	15.54	20.65	▁▆▇▃▁
Bioquimica_05	18	0.94	140.83	15.09	99.88	129.72	141.38	152.12	187.37	▁▆▇▃▁
Bioquimica_06	54	0.81	86.32	13.75	55.96	77.21	86.04	95.10	142.74	▃▇▅▁▁
Bioquimica_07	34	0.88	3.58	1.15	0.83	3.08	3.74	4.30	6.56	▂▂▇▃▁
Bioquimica_08	16	0.94	84.37	12.35	43.50	76.65	83.96	91.39	126.11	▁▃▇▃▁
Bioquimica_09	57	0.80	1.13	0.29	0.51	0.91	1.12	1.33	2.09	▃▇▆▂▁
Bioquimica_10	12	0.96	30.09	3.88	18.24	27.48	30.34	32.51	40.62	▁▃▇▆▁
Bioquimica_11	14	0.95	33.84	3.66	21.90	31.56	33.71	36.08	43.57	▁▃▇▆▂
Escala_01	32	0.89	103.18	122.93	0.00	21.23	54.69	137.77	892.17	▇▂▁▁▁
Escala_02	0	1.00	18.50	8.08	6.17	11.19	18.43	23.58	41.45	▇▇▇▃▁
Escala_03	0	1.00	18.65	8.24	5.79	12.32	17.76	24.00	46.96	▇▇▅▂▁
Genetica_01	160	0.44	0.80	0.60	0.00	0.29	0.64	1.23	2.99	▇▅▃▁▁
Genetica_02	170	0.40	0.78	0.55	0.01	0.33	0.73	1.05	2.27	▇▇▅▂▂
Genetica_03	178	0.38	0.76	0.56	0.00	0.32	0.67	1.09	2.22	▇▇▅▂▂
Genetica_04	168	0.41	0.84	0.55	0.01	0.42	0.70	1.20	2.60	▇▇▅▂▁
Genetica_05	145	0.49	20.19	214.06	0.00	0.84	1.82	3.17	2534.81	▇▁▁▁▁
Genetica_06	178	0.38	0.84	0.53	0.03	0.37	0.77	1.21	2.36	▇▇▆▂▁
Genetica_07	137	0.52	0.83	0.59	0.00	0.32	0.72	1.25	2.35	▇▆▅▃▂

#--- Valores perdidos
naniar::vis_miss(datos)

naniar::gg_miss_upset(datos, nsets = 10, nintersects = 50)

4 ¿De qué se compone tidymodels?

El flujo de trabajo para producir un modelo en tidymodels es el siguiente:

La idea detrás consiste en una sola sintaxis que resuma todo este flujo de trabajo.

Particion de datos
Pre-procesamos los datos de entrenamiento.
Buscamos los híper-parámetros ideales para nuestros datos en el modelo
Ajustamos el modelo.
Comprobamos con las métricas.
Vemos su fnciona con unos datos que el modelo no haya visto.

Como hemos visto en el ejemplo anterior, con cambiar una sola línea línea, podemos cambiar el tipo de modelo. Y con entender unos pocos puntos clave de todos los nexos, podemos llevar a cabo flujos de trabajo muy eficientes con muy poco coste de mantenimiento y poco cambio de sintaxis.

Para la creación de este “cosmos”, tidymodels utiliza las siguientes librerías:

rsample: Funciones para crear distintos tipos de remuestreos.
recipes: Preprocesador de datos, simplifica el proceso de preparado de datos para cualquier modelo.
parsnip: Interfaz para ajustar modelos, hasta 155 tipos a fecha de este html.
yarstick: Creación de métricas para la evaluación de un modelo.
workflows: Cohesión interna de tidymodels.
tune y dials: Creación y gestión de híper-parámetros de cualquier modelo presente en tidymodels.

5 Preparando la partición de datos con Rsampling

Para comenzar a crear un modelo de aprendizaje automático es necesario tener datos, fundamental. Estos datos, suelen ser difíciles de conseguir, requieren tiempo, esfuerzo y dinero. Su manejo debe ser optimizado con el fin de abaratar los costes de recopilarles y mejorar los resultados. Además tener en cuenta que el propósito de un modelo de aprendizaje automático en general es que se pueda usar para poder predecir que va a pasar con unos datos que no tenga, o que tenga en un futuro.

Aquí entre el punto del remuestreo. Crearemos un partición de la muestra en unos datos para “entrenar” el modelo y la otra para “evaluarlo”, para ver si tal y como ha quedado calibrado el modelo es útil o no.

Por supuesto, la calidad de esta partición estará ligada a la calidad de los datos. Si los datos recopilados no son suficientemente representativos, el modelo no predecirá bien cuando traigamos casos que no haya visto, aún incluso si la evaluación del modelo ha sido buena.

Básicamente queremos esto:

La librería dentro de tidymodels que nos dará los conjuntos de datos para hacer esto es rsample.

Algunas de las funciones importantes dentro de este conjunto:

initial_split(): Principal función dentro de la librería.Crea un objeto “split” sencillo para hacer una partición de datos

# 🔴 La función principal de Rsample, hace una partición de un conjunto de  🔴
initial_split(datos, prop = 0.6 )

## <Training/Testing/Total>
## <171/114/285>

bootstraps(): crea un tibble con n cantidad de objetos “split” y su correspondiente identificador. Estas muestras son de tipo bootstrap, es decir, muestras en las que un individuos puede aparecer más de una vez en el conjunto de datos, tanto en entrenamiento como en evaluación.

rsample::bootstraps(datos,5)

## # Bootstrap sampling 
## # A tibble: 5 × 2
##   splits            id        
##   <list>            <chr>     
## 1 <split [285/99]>  Bootstrap1
## 2 <split [285/103]> Bootstrap2
## 3 <split [285/113]> Bootstrap3
## 4 <split [285/104]> Bootstrap4
## 5 <split [285/113]> Bootstrap5

rolling_origin(): Crea un tibble con “n” cantidad de objetos “split” y su correspondiente identificador. En este caso, los remuestreos no son aleatorios y cada partición contienen datos que consecutivos. Es decir, se comienza por la fila 1 y se hacen bloqueen a cantidades iguales. Se usa sobretodo en la partición de datos de tiempo o cuando se está estudiando la creación de una partición optimizada entre el conjunto de entrenamiento y de evaluación.

rolling_origin(datos, initial = 20, assess = 20, skip = 40, cumulative = FALSE)

## # Rolling origin forecast resampling 
## # A tibble: 6 × 2
##   splits          id    
##   <list>          <chr> 
## 1 <split [20/20]> Slice1
## 2 <split [20/20]> Slice2
## 3 <split [20/20]> Slice3
## 4 <split [20/20]> Slice4
## 5 <split [20/20]> Slice5
## 6 <split [20/20]> Slice6

ggpubr::ggarrange(
  rolling_origin(datos, initial = 20, assess = 20, skip = 40, cumulative = FALSE) %>% 
    tidy() %>% 
    ggplot(aes(x = Resample, y = factor(Row), fill = Data)) +
    geom_tile() + 
  scale_fill_manual( values = c('cyan','orange'))
  ,
  rolling_origin(datos, initial = 50, assess = 20, skip = 70, cumulative = FALSE) %>% 
    tidy() %>% 
    ggplot(aes(x = Resample, y = factor(Row), fill = Data)) +
    geom_tile()  + 
  scale_fill_manual( values = c('cyan','orange'))
  ,
  rolling_origin(datos, initial = 20, assess = 20, skip = 40, cumulative = T) %>% 
    tidy() %>% 
    ggplot(aes(x = Resample, y = factor(Row), fill = Data)) +
    geom_tile()  + 
  scale_fill_manual( values = c('cyan','orange'))
  ,
  rolling_origin(datos, initial = 50, assess = 50, skip = 70, cumulative = T) %>% 
    tidy() %>% 
    ggplot(aes(x = Resample, y = factor(Row), fill = Data)) +
    geom_tile()  + 
  scale_fill_manual( values = c('cyan','orange'))
  ,
  common.legend = T, legend = 'bottom'
)

5.1 Estratificar

Una de las cosas que pasan habitual a la hora de crear una partición de datos es que, al ser aleatorias, uno de los conjuntos contenga más datos de un tiempo de otro tipo.

Idealmente querríamos que el conjunto de entrenamiento y que el conjunto de evaluación contuvieran la misma proporción de las categorías.

Para ello, en varias de las funciones de rsample se cuenta con el argumento strata que permite hacer esto de forma estratificada para una variable o variables de interés, es este caso, para nuestros datos y nuestra variable conjunta:

set.seed(4147)
# Splits ----

## Sin estratificar 
Split_datos_no_strat <- initial_split(datos, prop = 0.70 )
# Put 3/4 of the data into the training set 
# 🔴 Initial_split con datos estratificados 🔴
## Estratificados 
Split_datos_strat <- initial_split(datos, prop = 0.70, strata = Var_respuesta)

ggpubr::ggarrange(
  ggpubr::ggarrange(
    training(Split_datos_no_strat) %>% 
      group_by(Var_respuesta) %>% 
      summarise(recuento = n()) %>% 
      ungroup() %>% 
      mutate(prop_var_respuesta = recuento / sum(recuento) ) %>% 
      ggplot(aes(x = Var_respuesta, y = recuento, fill = Var_respuesta)) +
      geom_bar(stat = "identity") +
      geom_text(aes(label = scales::percent(prop_var_respuesta)), position = position_stack(vjust = 0.5), size = 4) +
    scale_fill_manual( values = c("azure4", "azure3"))
    ,
    testing(Split_datos_no_strat) %>% 
      group_by(Var_respuesta) %>% 
      summarise(recuento = n()) %>% 
      mutate(prop_var_respuesta = recuento / sum(recuento) ) %>% 
      ungroup() %>% 
      ggplot(aes(x = Var_respuesta, y = recuento, fill = Var_respuesta)) +
      geom_bar(stat = "identity") +
      geom_text(aes(label = scales::percent(prop_var_respuesta)), position = position_stack(vjust = 0.5), size = 4) +
      scale_fill_manual( values = c("azure4", "azure3")),
    common.legend = T, legend = 'bottom'
  ),
  ggpubr::ggarrange(
    training(Split_datos_strat) %>% 
      group_by(Var_respuesta) %>% 
      summarise(recuento = n()) %>% 
      ungroup() %>% 
      mutate(prop_var_respuesta = recuento / sum(recuento) ) %>% 
      ggplot(aes(x = Var_respuesta, y = recuento, fill = Var_respuesta)) +
      geom_bar(stat = "identity") +
      geom_text(aes(label = scales::percent(prop_var_respuesta)), position = position_stack(vjust = 0.5), size = 4) +
      scale_fill_manual( values = c("#E86EB7", "#108064"))
    ,
    testing(Split_datos_strat) %>% 
      group_by(Var_respuesta) %>% 
      summarise(recuento = n()) %>% 
      mutate(prop_var_respuesta = recuento / sum(recuento) ) %>% 
      ungroup() %>% 
      ggplot(aes(x = Var_respuesta, y = recuento, fill = Var_respuesta)) +
      geom_bar(stat = "identity") +
      geom_text(aes(label = scales::percent(prop_var_respuesta)), position = position_stack(vjust = 0.5), size = 4) +
      scale_fill_manual( values = c("#E86EB7", "#108064")), 
    common.legend = T,
    legend = 'bottom'
  ),
  labels = c('No estratificados', 'Estratificados'), label.x = 0.1
) %>%  
  plot()

Aclaración: los colores elegidos para el gráfico son de mi escena preferida de Spiderman, into the spiderverse. dejad siempre “huevos de Pascua” en los trabajos, hace que las cosas sean más emocionantes :P

5.2 Deshacer la partición

Para esto usaremos las funciones training() y test()

# 🔴 Deshacemos las particiones con las funciones training y testing 🔴

Training_datos <- training(Split_datos_strat)
Testing_datos  <- testing(Split_datos_strat)

knitr::kable(head(Training_datos,10))

Edad	Cognicion_01	Cognicion_02	Cognicion_03	Cognicion_04	Cognicion_05	Bioquimica_01	Bioquimica_02	Bioquimica_03	Bioquimica_04	Bioquimica_05	Bioquimica_06	Bioquimica_07	Bioquimica_08	Bioquimica_09	Bioquimica_10	Bioquimica_11	Escala_01	Escala_02	Escala_03	Genetica_01	Genetica_02	Genetica_03	Genetica_04	Genetica_05	Genetica_06	Genetica_07
14	66.78	NA	NA	NA	NA	NA	54.33	2.56	16.07	161.64	NA	NA	99.47	0.92	30.14	34.27	133.149284	18.116769	31.680105	1.6247538	NA	0.9775707	NA	0.9699578	NA	1.1983120
23	108.18	84.49	82.39	183.52	166.31	101.24	3.66	3.18	20.29	148.58	100.23	3.53	81.72	1.11	35.75	34.03	68.943373	12.183308	15.670177	NA	1.797438	NA	NA	0.1373531	1.1481697	NA
26	0.00	NA	NA	NA	NA	NA	129.41	5.99	14.98	145.04	94.53	4.20	78.76	0.87	32.57	26.89	2.898493	11.236520	6.628174	NA	NA	NA	1.6622978	2.6878823	1.6836672	0.1001570
14	52.88	NA	NA	NA	NA	NA	NA	2.06	13.98	132.94	105.75	NA	82.49	NA	34.73	29.83	99.600374	25.026820	12.949102	0.1247104	NA	NA	0.6977428	1.2486785	0.4195379	0.3760827
27	109.25	88.04	103.31	143.93	144.97	100.08	8.09	4.49	19.09	187.37	84.47	3.87	91.97	1.26	30.74	34.95	0.000000	16.798451	7.607729	0.8493442	NA	0.8781136	1.9962258	NA	NA	NA
30	105.60	104.09	86.45	149.99	105.02	93.65	7.20	4.51	15.65	145.00	90.13	3.57	90.67	1.05	28.98	31.50	51.087736	6.731747	10.894387	NA	1.349131	0.3128530	NA	1.9916249	1.2268437	NA
12	91.73	NA	NA	NA	NA	NA	30.29	3.15	13.94	134.92	103.40	3.97	96.92	1.27	21.42	29.26	160.752695	17.851526	18.519326	0.8739840	NA	0.5427362	NA	1.2463796	0.8887141	0.7917460
17	91.97	95.35	80.65	172.32	122.78	68.50	40.55	2.77	12.49	145.84	74.00	4.01	86.12	1.16	27.01	33.90	15.319511	14.058768	17.484487	0.2352628	NA	0.3039570	NA	NA	NA	0.9261596
14	69.65	NA	NA	NA	NA	NA	55.17	4.87	11.98	165.35	91.39	5.73	52.67	1.21	29.25	32.46	76.930164	22.999587	14.681546	NA	NA	NA	NA	NA	NA	NA
34	111.50	58.99	101.17	181.40	119.87	106.90	9.96	3.19	14.47	159.95	102.97	3.98	85.87	0.91	30.43	29.18	7.373436	8.195322	6.718927	1.2862183	NA	NA	NA	NA	NA	0.2123626

knitr::kable(head(Testing_datos,10))

Var_respuesta	Edad	Cognicion_01	Cognicion_02	Cognicion_03	Cognicion_04	Cognicion_05	Bioquimica_01	Bioquimica_02	Bioquimica_03	Bioquimica_04	Bioquimica_05	Bioquimica_06	Bioquimica_07	Bioquimica_08	Bioquimica_09	Bioquimica_10	Bioquimica_11	Escala_01	Escala_02	Escala_03	Genetica_01	Genetica_02	Genetica_03	Genetica_04	Genetica_05	Genetica_06	Genetica_07
0	14	101.82	NA	NA	NA	NA	NA	NA	5.65	13.65	125.71	NA	NA	84.29	1.48	30.35	35.68	0.00000	28.397013	7.363183	NA	0.9665678	NA	NA	1.292175	0.2417377	NA
0	17	80.41	85.53	55.68	124.52	159.70	62.37	NA	5.65	17.83	154.53	NA	NA	99.42	1.34	33.99	33.01	138.73119	16.769985	17.764831	0.4550571	NA	NA	NA	NA	NA	0.0044324
0	14	91.00	79.76	94.28	140.68	100.07	86.33	26.24	3.18	13.09	135.88	81.53	5.75	109.45	0.95	30.44	34.21	NA	21.382220	7.061147	NA	NA	NA	NA	NA	NA	NA
0	21	86.99	91.33	93.65	115.00	123.01	66.84	10.23	3.48	11.97	132.15	90.66	2.58	68.60	1.04	30.10	34.00	48.83394	7.785552	20.442332	NA	NA	NA	NA	NA	NA	NA
0	16	114.25	NA	NA	NA	NA	NA	8.14	5.68	12.25	149.89	NA	NA	91.22	1.28	35.86	40.11	NA	18.489371	10.365972	1.5581016	NA	NA	1.159770	NA	0.6242926	0.0579575
0	13	74.80	NA	NA	NA	NA	NA	93.32	4.77	8.88	164.33	NA	4.33	84.49	1.07	23.16	27.76	NA	18.749210	9.068513	NA	NA	NA	NA	NA	NA	NA
1	17	84.13	100.23	69.62	137.84	91.24	78.76	28.80	4.22	14.38	151.67	NA	5.23	91.19	1.75	18.24	37.76	113.35402	20.327740	21.101902	NA	NA	1.180101	NA	NA	1.0568568	NA
0	31	128.68	NA	99.23	173.48	144.55	101.62	21.34	5.61	14.16	147.96	104.95	4.06	91.39	1.23	31.89	32.58	118.10758	17.200513	14.212995	NA	NA	NA	0.972422	1.642294	NA	0.5527762
0	20	110.69	88.69	72.28	143.15	193.88	86.83	25.79	4.50	13.12	134.66	68.58	3.56	81.48	1.11	28.27	31.64	225.00452	17.789336	20.238253	NA	NA	NA	NA	NA	NA	NA
0	16	83.73	NA	NA	NA	NA	NA	23.32	4.16	12.12	137.10	83.00	5.53	71.56	0.99	32.46	35.09	133.40626	20.677244	15.434748	NA	NA	NA	NA	NA	NA	NA

5.3 folds, creando conjuntos de validación

# 🔴 Creamos Folds 🔴
Folds_training <- vfold_cv(Training_datos, v = 5 ,strata = Var_respuesta)

6 Pre-procesamiento de datos con `parsnip`

El pre-procesamiento es clave en el proceso de creación de un modelos de machine learning. Cuando queremos hacer un dataframe para modelizar en aprendizaje automático, no viene todo arreglado. Incluso con una recopilación de datos buena puede que necesitemos algunos arreglos. Algunas de las causas de esto pueden ser:

Hay modelos que no aceptan valores perdidos. Es necesario que imputemos (rellenar con valores perdidos).
Dependiendo de cómo queramos escalar el modelo o a qué pregunta estemos contestando, los datos deberán estar normalizados o no.
Re-equilibrar grupos disparejos aumenta la eficiencia de predicción, aunque hay que tener cuidado de cuando usarlo.
Reconvertir/transformar variables para maximizar el rendimiento del modelo.

Hacer todo estos procesos puede ser engorroso, pero podemos salir del paso fácilmente gracias por ejemplo de dplyr. Sin embargo ¿Cómo hacemos cuando queramos poner un modelo a grane scala y entren datos a cada momento? ¿Debemos ejecutar el pre-procesamiento una y otra vez?.

La clave para resolver esto está en el paquete recipes. Esta librería trae consigo todo un conjunto de funciones que permiten un pre-procesamiento fino, fino.

6.1 ¿Cómo cocinamos una receta?

Para poder elaborar un flujo completo de pre-procesamiento debemos seguir los siguuientes pasos:

recipe(): la base. Requiere de qué fórmula queremos para el modelo y de qué conjunto de datos partimos. Los datos aquí incorporados serán los datos de entrenamiento
step(): el proceso que queramos aplicar a nuestros datos, cualquier transformación de datos va aquí.
prep(): la receta se prepara, es decir se aplica al conjunto de entrenamiento y se abstrae para un modelado cualquiera
bake(): pone el pre-procesamiento en producción para un conjunto de datos cualquiera.

Hagamos un ejemplo básico: queremos una receta que sencillamente impute los datos faltantes por la media correspondiente a cada variable. Para eso vamos con ejemplo sencillo. La siguiente receta dejará los datos preparados para un modelo de aprendizaje automático, pero sin aplicarle ninguna transformación.

Training_datos %>% 
  # 🔴 1. Recipe (crea la formula generalizada del modelo) 
  recipe(  
    as.formula("Var_respuesta ~ Cognicion_01 +Cognicion_02+Escala_01+Escala_02"),
    data = . , 
    strata = Var_respuesta )  %>% 
  # 🔴 2. step
  # 🟡 En este caso no aplicamos ninguina trasnformación 
  # Step_log
  # 🔴 3. Prep (asienta la receta y generalizala) 
  prep() %>% 
  # 🔴 4. Bake (pon la receta en producción/comprueba que tal ha funcionado)
  bake(., new_data = NULL) %>% 
  # mostramos resultados
  head(.,10) %>% 
  knitr::kable()

Cognicion_01	Cognicion_02	Escala_01	Escala_02
66.78	NA	133.149284	18.116769
108.18	84.49	68.943373	12.183308
0.00	NA	2.898493	11.236520
52.88	NA	99.600374	25.026820
109.25	88.04	0.000000	16.798451
105.60	104.09	51.087736	6.731747
91.73	NA	160.752695	17.851526
91.97	95.35	15.319511	14.058768
69.65	NA	76.930164	22.999587
111.50	58.99	7.373436	8.195322

Training_datos %>% 
  # 🔴 1. Recipe (crea la formula generalizada del modelo) 
  recipe(  
    as.formula("Var_respuesta ~ Cognicion_01 +Cognicion_02+Escala_01+Escala_02"),
    data = . , 
    strata = Var_respuesta )  %>% 
  # 🔴 2. Step (haz la transofrmación que requieras) 
  # 🟡 Este step será el logaritmo de todas las variables predictoras númericas.
  step_log(all_numeric(), -all_outcomes()) %>% 
  # 🔴 3. Prep (asienta la receta y generalizala) 
  prep() %>% 
  # 🔴 4. Bake (pon la receta en producción/comprueba que tal ha funcionado)
  bake(., new_data = NULL) %>% 
  # mostramos resultados
  head(.,10) %>% 
  knitr::kable()

Cognicion_01	Cognicion_02	Escala_01	Escala_02
4.201404	NA	4.891471	2.896838
4.683796	4.436633	4.233286	2.500067
-Inf	NA	1.064191	2.419169
3.968025	NA	4.601166	3.219948
4.693639	4.477791	-Inf	2.821287
4.659658	4.645256	3.933544	1.906835
4.518849	NA	5.079867	2.882089
4.521462	4.557554	2.729127	2.643246
4.243483	NA	4.342898	3.135476
4.714025	4.077368	1.997884	2.103563

¡WALA! Datos pre-procesados y listos.

A partir de aquí en cielo es el límite ya que hay steps para todo, recomiendo leer la librería recipes y ver que tienen disponible. Al final de esta sección dejaremos lista una receta como ejemplo

6.2 Un par de ejemplos más

6.2.1 Imputación

Hacer una imputación simple suele no ser la mejor vía de hacerlo. lo ideal siempre es imputación múltiple. Es decir, una imputación que rellena los valores que faltan con números plausibles derivados de las distribuciones y relaciones entre las variables observadas en el conjunto de datos.

En tidymodels se puede hacer imputación múltiple, y valga la redundancia, de múltiples formas. Por ejemplo, podemos imputar una variable usando una regresión lineal a partir de otras. Hagamos una prueba con la variable Cognicion_02 a partir de Cognicion_01 y Escala_02

Training_datos %>% 
  # 🔴 1. Recipe (crea la formula generalizada del modelo) 
  recipe(  
    as.formula("Var_respuesta ~ Cognicion_01 +Cognicion_02+Escala_01+Escala_02"),
    data = . , 
    strata = Var_respuesta )  %>% 
  # 🔴 2. Step (haz la transofrmación que requieras) 
  # 🟡 Este step será la imputación por regresiópn lineal con otras covariables
  step_impute_linear( Cognicion_02, impute_with = imp_vars(Cognicion_01, Escala_02) ) %>%
  # 🔴 3. Prep (asienta la receta y generalizala) 
  prep() %>%
  # 🔴 4. Bake (pon la receta en producción/comprueba que tal ha funcionado)
  bake(., new_data = NULL) %>%   
  # mostramos resultados
  head(.,10) %>% 
  knitr::kable()

Cognicion_01	Cognicion_02	Escala_01	Escala_02
66.78	95.30378	133.149284	18.116769
108.18	84.49000	68.943373	12.183308
0.00	105.71896	2.898493	11.236520
52.88	97.79188	99.600374	25.026820
109.25	88.04000	0.000000	16.798451
105.60	104.09000	51.087736	6.731747
91.73	91.30365	160.752695	17.851526
91.97	95.35000	15.319511	14.058768
69.65	95.03224	76.930164	22.999587
111.50	58.99000	7.373436	8.195322

Esto va genial cuando conocemos la estructura de los datos y sabemos con certeza la dependencia de unas variables con otras para imputarlas fácilmente, en lugar de confiar ciegamente en algoritmos de imputación.

6.2.2 PCA

Como parte del pre-procesamiento, a veces puede ser necesario convertir ciertas variables a PCA para usarlas como componentes. No entraré en su uso, solo especificar la forma de hacerlo.

Training_datos %>% 
  # 🔴 1. Recipe (crea la formula generalizada del modelo) 
  recipe(  
    as.formula("Var_respuesta ~ Cognicion_01 + Cognicion_02 + Cognicion_03 + Escala_01 + Escala_02  + Escala_03"),
    data = . , 
    strata = Var_respuesta )  %>% 
  # 🔴 2. Step (haz la transofrmación que requieras) 
  # 🟡 Este step será la imputación por regresiópn lineal mcon otras covariables
  step_impute_bag(all_numeric()) %>% 
  step_pca( Cognicion_01 , Cognicion_02 , Cognicion_03, num_comp = 2,  id = "pca") %>%
  # 🔴 3. Prep (asienta la receta y generalizala) 
  prep() %>%
  # 🔴 4. Bake (pon la receta en producción/comprueba que tal ha funcionado)
  tidy(id = "pca") %>%
  # mostramos resultados
  head(.,10) %>% 
  knitr::kable()

terms	value	component	id
Cognicion_01	-0.5964794	PC1	pca
Cognicion_02	-0.6133048	PC1	pca
Cognicion_03	-0.5177543	PC1	pca
Cognicion_01	0.7421220	PC2	pca
Cognicion_02	-0.6671338	PC2	pca
Cognicion_03	-0.0647105	PC2	pca
Cognicion_01	0.3057241	PC3	pca
Cognicion_02	0.4228354	PC3	pca
Cognicion_03	-0.8530785	PC3	pca

6.2.3 Re-remuestrear

ADASYN, o Adaptive Synthetic Sampling, es un algoritmo bastante utilizado cuando tenemos conjuntos de datos desequilibrados.

Los conjuntos de datos desequilibrados se producen cuando la distribución de clases en los datos es desigual, lo que dificulta el entrenamiento de modelos que pueden estar sesgados hacia la clase mayoritaria.

https://www.researchgate.net/publication/224330873_ADASYN_Adaptive_Synthetic_Sampling_Approach_for_Imbalanced_Learning

# 🔴 Libreria que permite instalar algoritmos e sobremuestreo en tidymodels.
pacman::p_load(themis)

Training_datos_adasyn <- Training_datos %>% 
  # 1. Recipe (crea la formula generalizada del modelo) 
  recipe(  
    as.formula("Var_respuesta ~ Cognicion_01 +Cognicion_02+Escala_01+Escala_02"),
    data = . , 
    strata = Var_respuesta )  %>% 
  #  2. Step (haz la transofrmación que requieras) 
  #  step de  imputación por arboles manteniendo la estrucura de todas las covariables
  step_impute_bag(all_numeric()) %>% 
  # 🔴 step para sobremuestrear la el nivel menor en la variable respuesta.
  step_adasyn(Var_respuesta, over_ratio = 1, neighbors = 10) %>% 
  #  3. Prep (asienta la receta y generalizala) 
  prep() %>%
  #  4. Bake (pon la receta en producción/comprueba que tal ha funcionado)
  bake(., new_data = NULL)

# mostramos resultados
Training_datos_adasyn %>% 
  head(.,10) %>% 
  knitr::kable()

Cognicion_01	Cognicion_02	Escala_01	Escala_02
66.78	102.05505	133.149284	18.116769
108.18	84.49000	68.943373	12.183308
0.00	99.58442	2.898493	11.236520
52.88	97.03067	99.600374	25.026820
109.25	88.04000	0.000000	16.798451
105.60	104.09000	51.087736	6.731747
91.73	89.86937	160.752695	17.851526
91.97	95.35000	15.319511	14.058768
69.65	100.57019	76.930164	22.999587
111.50	58.99000	7.373436	8.195322

ggpubr::ggarrange(
  Training_datos %>% 
      group_by(Var_respuesta) %>% 
      summarise(recuento = n()) %>% 
      ungroup() %>% 
      mutate(prop_var_respuesta = recuento / sum(recuento) ) %>% 
      ggplot(aes(x = Var_respuesta, y = recuento, fill = Var_respuesta)) +
      geom_bar(stat = "identity") +
      geom_text(
        aes(label = scales::percent(prop_var_respuesta)), 
        position = position_stack(vjust = 0.5), size = 4) +
    scale_fill_manual( values = c("azure4", "azure3"))
    ,
    Training_datos_adasyn %>% 
      group_by(Var_respuesta) %>% 
      summarise(recuento = n()) %>% 
      mutate(prop_var_respuesta = recuento / sum(recuento) ) %>% 
      ungroup() %>% 
      ggplot(aes(x = Var_respuesta, y = recuento, fill = Var_respuesta)) +
      geom_bar(stat = "identity") +
      geom_text(
        aes(label = scales::percent(prop_var_respuesta)),
        position = position_stack(vjust = 0.5), size = 4) +
      scale_fill_manual( values = c("azure4", "azure3")),
    common.legend = T, legend = 'bottom'
  ,
  labels = c('Normal', 'Sobre Muestreo\n        Adasyn'), label.x = 0.1
) %>%  
  plot()

6.3 Caso práctico: receta

Por el momento dejaré una receta lista para continuar con el taller.

# Creamos la Fórmula para la receta a partir 
# de nombres de las variables y un poco de magia con paste
Receta_modelo_1_formula <- Training_datos %>% 
  dplyr::select(
   dplyr::matches('Cog'), dplyr::matches('escala')) %>% 
  names() %>% 
  paste(., collapse = ' + ') %>% 
  paste('Var_respuesta ~' ,.) %>% 
  as.formula() 

  # 🔴 1.) Iniciamos la receta
Receta_modelo_1 <- recipe(  
    formula = Receta_modelo_1_formula,
    data = Training_datos , 
    strata = Var_respuesta)  %>% 
  # 🔴 2.) Steps
    #  🟡 2.1) Eliminamos las variables que contengan más de un 20% de datos perdidos 
  step_filter_missing(all_predictors(),  threshold = 0.1) %>% 
    #  🟡 2.2) Imputamos las variables númericas con un algoritmo de bagged trees
  step_impute_bag(all_numeric(),-all_outcomes()) %>%
    # 2.3) Normalizamos los datos (restamos la media)
  step_normalize(all_numeric(),-all_outcomes() ) %>% 
    # 2.4) Escalamos los datos (reducimos a escala entre 0 y 1)
  step_scale(all_numeric(),-all_outcomes()) %>% 
    # 2.5) Convertimos a Dummy las variables factor (no hará nada por que nop tenemos variables factor)
  step_dummy(all_nominal(), -all_outcomes(), one_hot = T) %>% 
    # 2.6) Sobre muestreamos la variable para equilibrar grupos
  step_adasyn(Var_respuesta, over_ratio = 1, neighbors = 10) %>% 
  # 3.) Preparamos la receta
  prep()

Receta_modelo_1 %>% 
  bake(., new_data = NULL) %>% 
  head(., 10) %>% 
  knitr::kable()

Cognicion_01	Escala_02	Escala_03
-0.6426710	-0.0464912	1.5341196
0.7604689	-0.7638606	-0.3800329
-2.9059967	-0.8783295	-1.4610980
-1.1137735	0.7889501	-0.7053657
0.7967336	-0.2058789	-1.3439820
0.6730268	-1.4229671	-0.9510280
0.2029411	-0.0785598	-0.0393877
0.2110752	-0.5371131	-0.1631135
-0.5454002	0.5438529	-0.4982340
0.8729912	-1.2460174	-1.4502475

7 ¿Dónde está el modelo? Creando modelos con parnsip

Vale, dejamos atrás la parte del pre-procesamiento y empezamos la parte más machine leaninrg. Modelos. Modelos a la carrera.

La función del paquete parsnip dentro de tidymodels es propoconar una interfaz universal para diferentes métodos de modelado en R.

El ejemplo del principio del taller está basado en esta sección, parnsip es el encargado de que solo con cambiar unas pocas líneas puedas cambiar todo un modelo dentro de un flujo de aprendizaje automático.

A fecha de este taller, tidymodels cuenta con más de 155 modelos, incluyendo supervivencia y series de tiempo. Se pueden encontrar en: https://www.tidymodels.org/find/parsnip/

Supongamos que queremos empezar modelar, queremos usar una regresión logística, un clásico en este campo.Pero dentro de R hay varias librarías que usan la regresión logística, cada una con sus variaciones y reglas. Veamos que tenemos para esto:

# 🔴 Ver los posibles motores de un tipo de modelo
show_engines('logistic_reg')

## # A tibble: 8 × 2
##   engine    mode          
##   <chr>     <chr>         
## 1 glm       classification
## 2 glmnet    classification
## 3 LiblineaR classification
## 4 spark     classification
## 5 keras     classification
## 6 stan      classification
## 7 brulee    classification
## 8 h2o       classification

logistic_reg(
  engine = "glmnet",
  #Hyper-parametros           
  penalty = NULL, mixture = NULL, mode = 'classification')

## Logistic Regression Model Specification (classification)
## 
## Computational engine: glmnet

logistic_reg(
  engine = "glmnet",
  #Hyper-parametros           
  penalty = 0, mixture = 0, mode = 'classification')

## Logistic Regression Model Specification (classification)
## 
## Main Arguments:
##   penalty = 0
##   mixture = 0
## 
## Computational engine: glmnet

logistic_reg(
  engine = "glmnet",
  #Hyper-parametros           
  penalty = tune(), mixture = tune(), mode = 'classification')

## Logistic Regression Model Specification (classification)
## 
## Main Arguments:
##   penalty = tune()
##   mixture = tune()
## 
## Computational engine: glmnet

7.1 Caso práctico, modelo de Random Forest

Model_RandomForest <- 
  # 🔴 Especificamos el modelo que queremos en este caso un random forest
  rand_forest(
  # 🟡 Ajustamos los hiper-parametros 
    mtry = tune(),  trees = tune(),  min_n = tune()) %>% 
  # 🔴 Ponemos el motor, es decir la librería por la queremos que se ejecute el modelo
  set_engine("ranger", importance = "impurity") %>% 
  # 🔴 Ajustamos el modo del modelo, es decri si queremos regresion o clasifiacion 
  set_mode("classification")

8 Unificando todo: Workflow

El paquete de workflows permite combinar una receta y una especificación de modelo en un único objeto. Esto hace que todo el proceso de modelado, desde el pre-procesamiento de datos hasta el ajuste del modelo, esté en un mismo punto

WF_Receta_modelo_1_Random_forest <- 
  # 🔴 Activamos el workflow 
  workflow() %>% 
  # 🔴 Añadimos la receta
  add_recipe(Receta_modelo_1) %>% 
  # 🔴 Añadimos el modelo
  add_model(Model_RandomForest)

9 Eligiendo métricas con `yardstick`

El paquete yardstick, en concreto, se utiliza para las métricas de evaluación de modelos. Además

# 🔴 Seleccionamos las métricas que queremos para nuestro modelo
Modelo_Metricas <- metric_set(accuracy, j_index, precision, sensitivity, specificity, roc_auc, f_meas,recall, mcc)

También decir que es es posible crear tu propia métrica. Pero aún no he indagado en esta metodología. Se puede encontrar más en: https://yardstick.tidymodels.org/articles/custom-metrics.html

10 Eligiendo hiperparámetros: tune y dials

Los hiperparámetros son una pieza crucial en los modelos de aprendizaje automático. Estos son parámetros de configuración externos que no pueden aprenderse a partir de los datos, pero que influyen significativamente en el rendimiento del modelo. Un ejemplo sería el número de árboles de un random forest por ejemplo. A diferencia de los parámetros del modelo, que se aprenden durante el entrenamiento, los híperparámetros se establecen antes de que comience el entrenamiento de este.

La selección de los hiperparámetros adecuados es esencial para conseguir un modelo de aprendizaje automático generalizable y con buen rendimiento y es un proceso delicado. Los híper-parámetros desempeñan un papel crucial a la hora de controlar el equilibrio entre la sobreajuste y el infraajuste (overfitting y underfitting). El sobreajuste se produce cuando un modelo funciona bien con los datos de entrenamiento pero no consigue generalizarse a datos nuevos que no se han visto. Por otro lado, el infraajuste se produce cuando el modelo es demasiado simple y no puede captar los patrones subyacentes en los datos. Ajustar los hiperparámetros puede ayudar a encontrar el nivel adecuado de complejidad del modelo, haciendo un buen equilibrio entre la varianza del conjunto de entrenamiento y el de evaluación.

Los paquetes tune y dials son los que nos permiten hacer esta selección y control. Aquí se crea todo un subflujo de trabajo para el ajuste de híper-parámetros.

Creamos unas réplicas del conjunto de entrenamiento, para tener subconjunto de entrenamiento y de validación (que ya hicimos en el apartado de rsample con el objeto ““)
Elegimos un modelo con parsnip y dejamos sus híperparametros en abierto con tune().
Creamos una cuadrícula (grid) con diferentes posibles combinaciones de los híper-parámetros que dejamos en abierto en el paso 2.
Ponemos todo junto con la función tune_grid()

10.1 Montando un espacio de híper-parametros: familia grid_

grid_regular(): hará todas las combinaciones entre los rangos de la variable.
grid_random(): probará hiper parámetros aleatoriamente.
grid_max_entropy(): propondrá una combinacion de hiperpárametros que garanticen que se cubrirá todo el espectro, con fin de evitar mínimos locales en la conversión de algoritmos.
grid_latin_hypercube(): propondrá una matriz.

grid_regular() hará las combinaciones acorde al rango que le hayamos dado. Por defecto hace la combinación entre el mínimo, el medio y el máximo de cada valor.

grid_regular(
  mtry(range  = c(2, 10)),
  min_n(range = c(2, 10)),
  trees(range = c(2, 10)))

## # A tibble: 27 × 3
##     mtry min_n trees
##    <int> <int> <int>
##  1     2     2     2
##  2     6     2     2
##  3    10     2     2
##  4     2     6     2
##  5     6     6     2
##  6    10     6     2
##  7     2    10     2
##  8     6    10     2
##  9    10    10     2
## 10     2     2     6
## # ℹ 17 more rows

Podemos ampliar esto con el parámetro levels. Ahora será con el mínimo, el quantile 33, el cuantil 66 y el máximo.

grid_regular(
  mtry(range  = c(2, 10)),
  min_n(range = c(2, 10)),
  trees(range = c(2, 10)),
  levels = 4  )

## # A tibble: 64 × 3
##     mtry min_n trees
##    <int> <int> <int>
##  1     2     2     2
##  2     4     2     2
##  3     7     2     2
##  4    10     2     2
##  5     2     4     2
##  6     4     4     2
##  7     7     4     2
##  8    10     4     2
##  9     2     7     2
## 10     4     7     2
## # ℹ 54 more rows

grid_max_entropy(
  mtry(range  = c(1, 4)),
  min_n(range = c(10, 30)),
  trees(range = c(1, 1000)),
  size = 10
  )

## # A tibble: 10 × 3
##     mtry min_n trees
##    <int> <int> <int>
##  1     1    25    20
##  2     2    11   610
##  3     4    12     2
##  4     2    20   936
##  5     4    12   733
##  6     4    30   868
##  7     1    20   676
##  8     2    19   307
##  9     2    29   661
## 10     4    26   190

10.2 Ajustando los hiper-parametros: convergen `yardstick` con `dial` y `tune`

WF_hiper_parametros <- 
  # 🔴 Este proceso se hace con la función tune grid  🔴
  tune_grid(
  # Ponemos la receta nuestro modelo, que incuye el preprocesamiento y la fórmula
  object = WF_Receta_modelo_1_Random_forest,
  # 🟡 Ponemos los Folds del conjunto de entrenamiento,  
  # 🟡 Dentro del entrenamiento hará un sub entrenamiento y una sub validación
  resamples = Folds_training,
  #Ponemos un grid para que pruebe con diferentes híper-parámetros, ya con un rango pre-definido
  grid = grid_max_entropy(
    mtry(range  = c(1, 4)),
    min_n(range = c(10, 30)),
    trees(range = c(1, 1000)),
    size = 10),
  # Ponemos las métricas que hemos definido anteriormente
  metrics = Modelo_Metricas, 
  # Con esta opción podemos guardar las predicciones
  control = control_grid( save_pred = T)
)

WF_hiperparametros_coleccion <- 
  WF_hiper_parametros %>% 
  # 🔴 Esta es la funcion que retorna las métricas, importante tenerla en cuenta
  collect_metrics() 

WF_hiperparametros_coleccion %>% 
  head(.,20) %>% 
  knitr::kable()

mtry	trees	min_n	.metric	.estimator	mean	n	std_err	.config
2	886	10	accuracy	binary	0.6610256	5	0.0496335	Preprocessor1_Model01
2	886	10	f_meas	binary	0.7673926	5	0.0365573	Preprocessor1_Model01
2	886	10	j_index	binary	0.1832661	5	0.1253188	Preprocessor1_Model01
2	886	10	mcc	binary	0.1626371	5	0.1086564	Preprocessor1_Model01
2	886	10	precision	binary	0.8400397	5	0.0310967	Preprocessor1_Model01
2	886	10	recall	binary	0.7082661	5	0.0439273	Preprocessor1_Model01
2	886	10	roc_auc	binary	0.7213836	5	0.0469202	Preprocessor1_Model01
2	886	10	sensitivity	binary	0.7082661	5	0.0439273	Preprocessor1_Model01
2	886	10	specificity	binary	0.4750000	5	0.0918559	Preprocessor1_Model01
4	123	28	accuracy	binary	0.6760256	5	0.0501729	Preprocessor1_Model02
4	123	28	f_meas	binary	0.7747366	5	0.0350463	Preprocessor1_Model02
4	123	28	j_index	binary	0.2768145	5	0.1524681	Preprocessor1_Model02
4	123	28	mcc	binary	0.2346557	5	0.1298836	Preprocessor1_Model02
4	123	28	precision	binary	0.8684292	5	0.0392609	Preprocessor1_Model02
4	123	28	recall	binary	0.7018145	5	0.0393757	Preprocessor1_Model02
4	123	28	roc_auc	binary	0.7346900	5	0.0462352	Preprocessor1_Model02
4	123	28	sensitivity	binary	0.7018145	5	0.0393757	Preprocessor1_Model02
4	123	28	specificity	binary	0.5750000	5	0.1286954	Preprocessor1_Model02
2	81	16	accuracy	binary	0.6762821	5	0.0479575	Preprocessor1_Model03
2	81	16	f_meas	binary	0.7791279	5	0.0342760	Preprocessor1_Model03

autoplot(WF_hiper_parametros) +
  theme_bw()

WF_hiper_parametros %>% 
  collect_metrics() %>% 
  pivot_longer(cols = c('mtry','trees', 'min_n'), names_to = 'tipo_parametro', values_to = 'valor_parametro') %>% 
  ggplot(aes(valor_parametro,mean, color = mean)) +
  geom_point() +
  facet_grid(.metric~tipo_parametro, scales = 'free' ) +
  scale_color_viridis_b() +
  theme_bw()

WF_hiperparametros_mejor <- WF_hiper_parametros %>%  show_best("roc_auc",1)

WF_hiperparametros_mejor %>% 
  knitr::kable()

mtry	trees	min_n	.metric	.estimator	mean	n	std_err	.config
2	240	29	roc_auc	binary	0.7387727	5	0.0384544	Preprocessor1_Model06

11 La conversión de todo, `fit` y `predict`

Ya hemos hecho todos los procesos previos. Todo lo que a un modelo de machine learning le puede hacer falta. Finalmente, una vez que lo tenemos todo hilado, hay que proceder a entrenar el modelo y dejarlo listo.

Modelo_fitted <- WF_Receta_modelo_1_Random_forest %>% 
  # 🔴 Finalizamos el flujo de trabajo con la mejor combinación de Hiper paraemtros encontrada
  finalize_workflow(WF_hiperparametros_mejor) %>%
  # 🔴 finalmente ajustamos con fit
  fit(data = Training_datos)

Con fit ya por fin tenemos nuestro modelo, ya podemos ponernos a hacer lo que hace el aprendizaje automático: predecir. vamos a crear una observación nueva y a ver como la trabaja.

Nuevo_caso <- tibble(
  'Cognicion_01' = 123,
  'Cognicion_02' = 102,
  'Cognicion_03' = 89,
  'Cognicion_04' = 100,
  'Cognicion_05' = 179,
  'Escala_01'    = 12,
  'Escala_02'    = 22,
  'Escala_03'    = 20
)

tibble(
  # 🔴 La función predict de la librería stats tiene sinergias con tidymodels, neecsitra de estos argumentos:
  # - 🟡 El modelo ya ajustado (model_fitted) 
  # - 🟡 Un tibble/dataframe con la nueva prediccion (debe contener las mismas variables)
  # - 🟡 Un tipo de prediccion en los casos de clasificación, si queremos las probabiliadades de cada categoria o la categoría predicha  directamente
  predict(Modelo_fitted, Nuevo_caso, type = "prob" ),
  predict(Modelo_fitted, Nuevo_caso, type = "class")) %>% 
  knitr::kable()

.pred_0	.pred_1	.pred_class
0.9541702	0.0458298	0

tibble(predict(Modelo_fitted, Testing_datos, type = "prob" ), predict(Modelo_fitted, Testing_datos, type = "class")) %>% 
  head(.,10) %>% 
  knitr::kable()

.pred_0	.pred_1	.pred_class
0.7868158	0.2131842	0
0.5619997	0.4380003	0
0.8653084	0.1346916	0
0.4211667	0.5788333	1
0.8627070	0.1372930	0
0.9247541	0.0752459	0
0.0893169	0.9106831	1
0.9033712	0.0966288	0
0.9443273	0.0556727	0
0.7467458	0.2532542	0

11.1 `last_fit`

Crear las predicciones y todo es un sistema engorroso. Además de que se pierde mucha informaciíon por el camino, ya que todo sigue estando en objetos dispersos. Para ello los desarrolladores crearon last_fit.

La función last_fit permite tener todo unificado en un solo tibble, al tener que ajustarse con un objeto split, esta recibe directamente los datos de entrenamiento y los datos de evaluación (recibe datos de training y testing) y entonces crea sus métricas y sus predicciones.

Con last_fit() contamos con una serie de funciones para extraer todo lo que necesitemos de ellas:

extract_fit_engine() : permite extraer el modelo crudo como si no hubiera sido ajustado con tidymodels, de modo que si es un modelo explicativo puede ser tratado con normalidad
collect_metrics() : extraer un tibble con las métricas del modelo (puestas con metric set)
collect_predictions() : extrae un tibble con las predicciones de los datos de evaluación. (Por defecto el umbral de rendimiento en esta función es 0.5, y estoy buscando alguna forma de cambiar esto, pero se puede manipular el mismo tibble como cualquier otro, de modo que no es un drama.)
extract_workflow() : permite extraer todo el flujo del modelo, de incluyendo la receta de modo que podemos hacer predicciones:

Modelo_1_final <-  WF_Receta_modelo_1_Random_forest %>% 
  # 🔴 Finalizamos el flujo de trabajo con la mejor combinación de Hiper paraemtros encontrada
  finalize_workflow(WF_hiperparametros_mejor) %>%
  # 🔴 Ajustamos con last_fit, poniendo el split de datos inicial y las metricas para las predicciones
  last_fit(Split_datos_strat, metrics = Modelo_Metricas)

Modelo_1_final

## # Resampling results
## # Manual resampling 
## # A tibble: 1 × 6
##   splits           id               .metrics .notes   .predictions .workflow 
##   <list>           <chr>            <list>   <list>   <list>       <list>    
## 1 <split [198/87]> train/test split <tibble> <tibble> <tibble>     <workflow>

Modelo_1_final %>%
  # 🔴 dentro del flujo de trabajo están la receta y el modelo ya ajustado
  extract_workflow() %>% 
  # de modo que lo podemos usar para hacer predicciones
  predict(Testing_datos, type = "prob" )

## # A tibble: 87 × 2
##    .pred_0 .pred_1
##      <dbl>   <dbl>
##  1  0.787   0.213 
##  2  0.562   0.438 
##  3  0.865   0.135 
##  4  0.421   0.579 
##  5  0.863   0.137 
##  6  0.925   0.0752
##  7  0.0893  0.911 
##  8  0.903   0.0966
##  9  0.944   0.0557
## 10  0.747   0.253 
## # ℹ 77 more rows

A partir de aquí podemos manipular el objeto fit de múltiples maneras múltiples.

Modelo_1_final %>%
  # 🔴 Con esta funcion podemos recopilar todas las metricas de un modelo
  collect_metrics() %>% 
  select(-.estimator, -.config) %>%  
  mutate(.estimate = round(.estimate,2)) %>% 
  knitr::kable()

.metric	.estimate
accuracy	0.71
j_index	0.39
precision	0.89
sensitivity	0.72
specificity	0.67
f_meas	0.80
recall	0.72
mcc	0.33
roc_auc	0.74

options(yardstick.event_first = FALSE)

Modelo_1_final %>%  collect_predictions() %>% 
  # 🔴 Esta funcion permite calcular rapidamente la curva roc de un modelo ya ajustado
  roc_curve(truth = Var_respuesta,  .pred_0 ) %>% 
  ggplot(aes( x = 1 - specificity ,y = sensitivity ) ) +
  geom_path() +
  geom_abline(intercept = 0 , slope = 1, linetype = 3) +
  labs(title = 'Cruva roc del modelo_1_final') +
  theme_bw()

Advertencia: en la función de curva roc, ved que puesto predicción del 0 en vez del 1. Esto es un error interno en la clasificación dentro del del paquete yardstick. Hay que tener cuidado sobre si estas funciones utilizan el primer / segundo / x nivel de su factor como el “evento”. Por defecto, yardstick elige el primer nivel de verdad como “evento” al calcular la curva roc. Podéis encontrar más de este suceso en el github de los desarrolladores: https://github.com/tidymodels/yardstick/issues/94

12 Rendimiento del umbral (threshold performance en modelos de clasificación)

El concepto de rendimiento del umbral (threshold performance) es crucial a la hora de evaluar y ajustar modelos de clasificación binaria.

El umbral es el valor de probabilidad por encima del cual un resultado predicho se considera de clase positiva. Por defecto, muchos modelos utilizan un umbral de 0,5. Sin embargo, en casi todos los casos, el umbral óptimo no estará otro punto. Este “óptimo lo tendremos que definir que en función de los objetivos y requisitos específicos del problema.

Ajustar el umbral permite elegir entre precisión y la recall. Bajar el umbral aumenta la sensibilidad (recall) pero disminuye la precisión, y viceversa. Básicamente, el umbral óptimo depende de la importancia relativa de los falsos positivos y los falsos negativos en una aplicación concreta.

En tidymodels está la librería probably con la función threshold_perf(), la cual nos permite, dadas las predicciones del modelo ver como se calcula este umbral.

pacman::p_load(probably)


Modelo_1_final_Threshold_performance <- 
  # Llamamos al modelo ya finallizado con "last_fit"
  Modelo_1_final %>% 
  # Coleccionamos sus predicciones
  collect_predictions() %>% 
  # 🔴 Ejecutamos el redimeinto de umbral con la siguiente función
  threshold_perf(
    # Le decimos que variable es la verdadera respuesta
    truth = Var_respuesta,  
    # le decimos que variable es la predicción (aviso: debido a un bugg en la version )
    .pred_1, 
    # ajustamos un rango en el que se evaluará el umbral en cada punto
    thresholds = seq(0.5,1,0.01) )

# displayment
Modelo_1_final_Threshold_performance %>% 
  head(., 10) %>% 
  mutate(.estimate = round(.estimate, 2)) %>% 
  knitr::kable()

.threshold	.metric	.estimator	.estimate
0.50	sensitivity	binary	0.28
0.51	sensitivity	binary	0.28
0.52	sensitivity	binary	0.28
0.53	sensitivity	binary	0.28
0.54	sensitivity	binary	0.28
0.55	sensitivity	binary	0.28
0.56	sensitivity	binary	0.26
0.57	sensitivity	binary	0.25
0.58	sensitivity	binary	0.23
0.59	sensitivity	binary	0.23

Modelo_1_final_Threshold_performance %>% 
  ggplot(aes(x = .threshold, y = .estimate, color = .metric)) +
  geom_line(size = 1.1) +
  scale_color_manual(values = c('purple','blue','pink')) +
  geom_vline(xintercept = 
    Modelo_1_final_Threshold_performance %>%  
      filter(.metric == 'j_index') %>%
      arrange(desc(.estimate)) %>% 
      slice(1) %>% 
      pull(.threshold), linetype = 'dashed', size = 1.5)

13 El cielo es el límite, `workflow_set()`

Bien, ya hemos aprendido a crear todo un flujo entero de machine learning con tidymodels, desde la base hasta su final. Muy útil todo. pero calibrar modelos de 1 en 1 es un auténtico petardazo. No me acaba de ser útil este taller. ¿Cómo me va a ser útil si tengo que probar unas 20 fórmulas con datos?

Lo entiendo perfectamente, a mi también me pasó. Los médicos me preguntaron que era mejor para tener en cuenta, si las escalas, las cogniciones, la bioquímica o la genética a la hora de diagnosticar al paciente.

Comprobémoslo…

13.1 Todas las recetas

Receta_modelo_Escalas_formula <- Training_datos %>% 
  dplyr::select(dplyr::matches('Esc')) %>% 
  names() %>% 
  paste(., collapse = ' + ') %>% 
  paste('Var_respuesta ~' ,.) %>% 
  as.formula() 

  # 🔴 1.) Iniciamos la receta
Receta_modelo_Escalas <- recipe(  
    formula = Receta_modelo_Escalas_formula,
    data = Training_datos , 
    strata = Var_respuesta)  %>% 
  # 🔴 2.) Steps
    #  🟡 2.1) Eliminamos las variables que contengan más de un 25% de datos perdidos 
  step_filter_missing(all_predictors(),  threshold = 0.25) %>% 
    #  🟡 2.2) Imputamos las variables númericas con un algoritmo de bagged trees
  step_impute_bag(all_numeric(),-all_outcomes()) %>%
    # 2.3) Normalizamos los datos (restamos la media)
  step_normalize(all_numeric(),-all_outcomes() ) %>% 
    # 2.4) Escalamos los datos (reducimos a escala entre 0 y 1)
  step_scale(all_numeric(),-all_outcomes()) %>% 
    # 2.5) Convertimos a Dummy las variables factor (no hará nada por que nop tenemos variables factor)
  step_dummy(all_nominal(), -all_outcomes(), one_hot = T) %>% 
    # 2.6) Sobre muestreamos la variable para equilibrar grupos
  step_adasyn(Var_respuesta, over_ratio = 1, neighbors = 10) %>% 
  # 3.) Preparamos la receta
  prep()



Receta_modelo_Cognicion_formula <- Training_datos %>% 
  dplyr::select(
   dplyr::matches('Cog')) %>% 
  names() %>% 
  paste(., collapse = ' + ') %>% 
  paste('Var_respuesta ~' ,.) %>% 
  as.formula() 

  # 🔴 1.) Iniciamos la receta
Receta_modelo_Cognicion <- recipe(  
    formula = Receta_modelo_Cognicion_formula,
    data = Training_datos , 
    strata = Var_respuesta)  %>% 
  # 🔴 2.) Steps
    #  🟡 2.1) Eliminamos las variables que contengan más de un 25% de datos perdidos 
  step_filter_missing(all_predictors(),  threshold = 0.25) %>%
    #  🟡 2.2) Imputamos las variables númericas con un algoritmo de bagged trees
  step_impute_mean(all_numeric(),-all_outcomes()) %>%
    # 2.3) Normalizamos los datos (restamos la media)
  step_normalize(all_numeric(),-all_outcomes() ) %>% 
    # 2.4) Escalamos los datos (reducimos a escala entre 0 y 1)
  step_scale(all_numeric(),-all_outcomes()) %>% 
    # 2.5) Convertimos a Dummy las variables factor (no hará nada por que nop tenemos variables factor)
  step_dummy(all_nominal(), -all_outcomes(), one_hot = T) %>% 
    # 2.6) Sobre muestreamos la variable para equilibrar grupos
  step_adasyn(Var_respuesta, over_ratio = 1, neighbors = 10) %>% 
  # 3.) Preparamos la receta
  prep()

bake(Receta_modelo_Cognicion, new_data = NULL)

## # A tibble: 316 × 4
##    Cognicion_01 Cognicion_03 Cognicion_04 Var_respuesta
##           <dbl>        <dbl>        <dbl> <fct>        
##  1       -0.643    -2.70e-15     2.51e-15 0            
##  2        0.760     3.13e- 1     1.13e+ 0 0            
##  3       -2.91     -2.70e-15     2.51e-15 0            
##  4       -1.11     -2.70e-15     2.51e-15 0            
##  5        0.797     1.64e+ 0     2.59e- 1 0            
##  6        0.673     5.70e- 1     3.93e- 1 0            
##  7        0.203    -2.70e-15     2.51e-15 0            
##  8        0.211     2.03e- 1     8.86e- 1 0            
##  9       -0.545    -2.70e-15     2.51e-15 0            
## 10        0.873     1.50e+ 0     1.09e+ 0 0            
## # ℹ 306 more rows

Receta_modelo_Bioquimicas_formula <- Training_datos %>% 
  dplyr::select(
   dplyr::matches('Bio')) %>% 
  names() %>% 
  paste(., collapse = ' + ') %>% 
  paste('Var_respuesta ~' ,.) %>% 
  as.formula() 

  # 🔴 1.) Iniciamos la receta
Receta_modelo_Bioquimica <- recipe(  
    formula = Receta_modelo_Bioquimicas_formula,
    data = Training_datos , 
    strata = Var_respuesta)  %>% 
  # 🔴 2.) Steps
    #  🟡 2.1) Eliminamos las variables que contengan más de un 25% de datos perdidos 
  step_filter_missing(all_predictors(),  threshold = 0.25) %>% 
    #  🟡 2.2) Imputamos las variables númericas con un algoritmo de bagged trees
  step_impute_mean(all_numeric()) %>%
    # 2.3) Normalizamos los datos (restamos la media)
  step_normalize(all_numeric(),-all_outcomes() )  %>% 
    # 2.4) Escalamos los datos (reducimos a escala entre 0 y 1)
  step_scale(all_numeric(),-all_outcomes()) %>% 
    # 2.5) Convertimos a Dummy las variables factor (no hará nada por que nop tenemos variables factor)
  step_dummy(all_nominal(), -all_outcomes(), one_hot = T) %>% 
    # 2.6) Sobre muestreamos la variable para equilibrar grupos
  step_adasyn(Var_respuesta, over_ratio = 1, neighbors = 10) %>% 
  # 3.) Preparamos la receta
  prep()

Receta_modelo_Popurri <- Training_datos %>%
  dplyr::select(
    Cognicion_01, Cognicion_02, Escala_01,Escala_02, Bioquimica_01,Bioquimica_02, Genetica_01,Genetica_02) %>% 
  names() %>% 
  paste(., collapse = ' + ') %>% 
  paste('Var_respuesta ~' ,.) %>% 
  as.formula() 

  # 🔴 1.) Iniciamos la receta
Receta_modelo_Popurri <- recipe(  
    formula = Receta_modelo_Popurri,
    data = Training_datos , 
    strata = Var_respuesta)  %>% 
  # 🔴 2.) Steps
    #  🟡 2.1) Eliminamos las variables que contengan más de un 5% de datos perdidos 
  step_filter_missing(all_predictors(),  threshold = 0.25) %>% 
    #  🟡 2.2) Imputamos las variables númericas con un algoritmo de bagged trees
  step_impute_bag(all_numeric(),-all_outcomes()) %>%
    # 2.3) Normalizamos los datos (restamos la media)
  step_normalize(all_numeric(),-all_outcomes() ) %>% 
    # 2.4) Escalamos los datos (reducimos a escala entre 0 y 1)
  step_scale(all_numeric(),-all_outcomes()) %>% 
    # 2.5) Convertimos a Dummy las variables factor (no hará nada por que nop tenemos variables factor)
  step_dummy(all_nominal(), -all_outcomes(), one_hot = T) %>% 
    # 2.6) Sobre muestreamos la variable para equilibrar grupos
  step_adasyn(Var_respuesta, over_ratio = 1, neighbors = 10) %>% 
  # 3.) Preparamos la receta
  prep()

13.2 `workflow_set()` para unificarlas todas

Ahora, esto son muchas recetas y mucho texto. Será difícil combinarlo todo, ¿no?…¡Qué va!. workflow_set() ya está pensado para esto. Funciona como el crossing() de la librería purrr o como el expand.grid() de la librería base, crear workflows con todas las combinaciones.

Receta_list <- list(
  'Modelo_Escalas'    = Receta_modelo_Escalas,
  'Modelo_Cognicion'  = Receta_modelo_Cognicion,
  'Modelo_Bioquimica' = Receta_modelo_Bioquimica,
  'Modelo_Popurri'    = Receta_modelo_Popurri
  )
  
Model_list <- list(
  RandomForest     = Model_RandomForest
  # ,XGBoost          = Model_XGBoost
  # ,Bagged_Trees     = Model_Bagged_Tree
  )

Wokflows_set <- workflow_set(
  preproc = Receta_list, 
  models = Model_list, 
  # 🔴 La opción Cross es para que haga todos los cruces de modelos con recetas
  cross = T)


Wokflows_set

## # A workflow set/tibble: 4 × 4
##   wflow_id                       info             option    result    
##   <chr>                          <list>           <list>    <list>    
## 1 Modelo_Escalas_RandomForest    <tibble [1 × 4]> <opts[0]> <list [0]>
## 2 Modelo_Cognicion_RandomForest  <tibble [1 × 4]> <opts[0]> <list [0]>
## 3 Modelo_Bioquimica_RandomForest <tibble [1 × 4]> <opts[0]> <list [0]>
## 4 Modelo_Popurri_RandomForest    <tibble [1 × 4]> <opts[0]> <list [0]>

Una vez hecho el workflow_set, se pone a trabajar con workflow_map.

Wokflows_set_map <- 
  Wokflows_set %>% 
  # 🔴 workflow maps es una funcion que permite ajustar múltiples flujos de trabajo
    workflow_map(
    resamples = Folds_training, 
    fn = "tune_grid",
    grid = grid_max_entropy(
      mtry(range  = c(1, 4)),
      min_n(range = c(10, 30)),
      trees(range = c(1, 1000)),
      size = 10),
    # verbose = TRUE, 
    metrics = Modelo_Metricas, 
    control = control_grid( save_pred = T),
    # 🔴 para garantizare replicabildiad, podemos fijar la semilla en el proceso
    seed = 2465)

Wokflows_set_map

## # A workflow set/tibble: 4 × 4
##   wflow_id                       info             option    result   
##   <chr>                          <list>           <list>    <list>   
## 1 Modelo_Escalas_RandomForest    <tibble [1 × 4]> <opts[4]> <tune[+]>
## 2 Modelo_Cognicion_RandomForest  <tibble [1 × 4]> <opts[4]> <tune[+]>
## 3 Modelo_Bioquimica_RandomForest <tibble [1 × 4]> <opts[4]> <tune[+]>
## 4 Modelo_Popurri_RandomForest    <tibble [1 × 4]> <opts[4]> <tune[+]>

AVISO, workflow_map puede paralelizarse, pero aún estoy descubriendo como integrarlo de forma efectiva.

Finalmente, podemos calibrar cada modelo individualmente y tener todo finiquitado. Ese proceso sería más eficiente con una función, desde luego. Pero se me acaba el tiempo para dedicar al taller :(

# Modelo Escalas ----
Modelo_Escalas_RandomForest_hyperparametros <- Wokflows_set_map %>% 
  # 🔴 Extraemos el flujo correspondiente al modelo que queremos 
  extract_workflow_set_result('Modelo_Escalas_RandomForest') %>%
  # 🔴 podemos ver una tabla con los mejores resultados, acorde a una metrica
  show_best(metric = 'sensitivity', n = 50) %>% 
  filter(row_number() == 1) %>% 
  select(mtry,trees,min_n,.config )

Modelo_Escalas_RandomForest_final <- Wokflows_set_map %>% 
  # 🔴 Extraemos el flujo correspondiente al modelo que queremos 
  extract_workflow('Modelo_Escalas_RandomForest') %>% 
  finalize_workflow(Modelo_Escalas_RandomForest_hyperparametros) %>%
  last_fit(Split_datos_strat, metrics = Modelo_Metricas )

# Modelo_Cognicion ----

Modelo_Cognicion_RandomForest_hyperparametros <- Wokflows_set_map %>% 
  extract_workflow_set_result('Modelo_Cognicion_RandomForest') %>%
  show_best(metric = 'sensitivity', n = 50) %>% 
  filter(row_number() == 1) %>% 
  select(mtry, trees, min_n, .config )

Modelo_Cognicion_RandomForest_final <- Wokflows_set_map %>% 
  extract_workflow('Modelo_Cognicion_RandomForest') %>% 
  finalize_workflow(Modelo_Cognicion_RandomForest_hyperparametros) %>%
  last_fit(Split_datos_strat, metrics = Modelo_Metricas )

# Modelo_Bioquimico ----

Modelo_Bioquimica_RandomForest_hyperparametros <- Wokflows_set_map %>% 
  extract_workflow_set_result('Modelo_Bioquimica_RandomForest') %>%
  show_best(metric = 'sensitivity', n = 50) %>% 
  filter(row_number() == 1) %>% 
  select(mtry, trees, min_n, .config )

Modelo_Bioquimica_RandomForest_final <- Wokflows_set_map %>% 
  extract_workflow('Modelo_Bioquimica_RandomForest') %>% 
  finalize_workflow(Modelo_Bioquimica_RandomForest_hyperparametros) %>%
  last_fit(Split_datos_strat, metrics = Modelo_Metricas )

# Modelo Popurri ----
Modelo_Popurri_RandomForest_hyperparametros <- Wokflows_set_map %>% 
  extract_workflow_set_result('Modelo_Popurri_RandomForest') %>%
  show_best(metric = 'sensitivity', n = 50) %>% 
  filter(row_number() == 1) %>% 
  select(mtry, trees, min_n, .config )

Modelo_Popurri_RandomForest_final <- Wokflows_set_map %>% 
  extract_workflow('Modelo_Popurri_RandomForest') %>% 
  finalize_workflow(Modelo_Popurri_RandomForest_hyperparametros) %>%
  last_fit(Split_datos_strat, metrics = Modelo_Metricas )

Una vez acabado todo, podemos ver todo lo necesario usando unos pocos verbos como collect() y las lógicas del tidyverse.

13.3 Métricas de todos los modelos

map2(
  list(
    Modelo_Escalas_RandomForest_final,
    Modelo_Cognicion_RandomForest_final,
    Modelo_Bioquimica_RandomForest_final,
    Modelo_Popurri_RandomForest_final),
  list('Modelo_Escala','Modelo_Cognicion','Modelo_Bioquimica','Modelo_Popurri'),
  ~ ..1 %>% 
    collect_metrics() %>% 
    mutate('Modelo'= ..2) %>% 
    select(.metric,.estimate,Modelo )) %>% 
  bind_rows() %>% 
  pivot_wider(names_from = Modelo, values_from = .estimate) %>% 
  knitr::kable()

.metric	Modelo_Escala	Modelo_Cognicion	Modelo_Bioquimica	Modelo_Popurri
accuracy	0.6436782	0.6091954	0.7471264	0.6666667
j_index	0.2222222	-0.0265700	0.0241546	0.2101449
precision	0.8518519	0.7868852	0.7974684	0.8448276
sensitivity	0.6666667	0.6956522	0.9130435	0.7101449
specificity	0.5555556	0.2777778	0.1111111	0.5000000
f_meas	0.7479675	0.7384615	0.8513514	0.7716535
recall	0.6666667	0.6956522	0.9130435	0.7101449
mcc	0.1855216	-0.0235126	0.0338612	0.1805788
roc_auc	0.7504026	0.5024155	0.5048309	0.6231884

13.4 Cruvas ROC todos los modelos

ROC_CURVE_training_todos_los_posibles_modelos <- Wokflows_set_map %>%  
  collect_predictions() %>%
  group_split(wflow_id) %>% 
  set_names(list('Escala','Cognicion','Bioquimica','Popurri')) %>% 
  map(~ .x %>% roc_curve(truth = Var_respuesta , .pred_0)) %>% 
  map2(., list('Escala','Cognicion','Bioquimica','Popurri'),
               ~ .x %>%  mutate(model = .y)) %>% 
         bind_rows() 

Plot_ROC_CURVE_training_todos_los_posibles_modelos <- ROC_CURVE_training_todos_los_posibles_modelos %>% 
  ggplot(aes( x = 1 - specificity,y = sensitivity, color = model ) ) +
  geom_path() +
  geom_abline(intercept = 0, slope = 1, linetype = 3) +
  theme_bw()
  


ROC_CURVE_training_modelos_finales <- list(
  list('Esc' ,Modelo_Escalas_RandomForest_hyperparametros$.config),
  list('Cog' ,Modelo_Cognicion_RandomForest_hyperparametros$.config),
  list('Bio' ,Modelo_Bioquimica_RandomForest_hyperparametros$.config),
  list('Pop' ,Modelo_Popurri_RandomForest_hyperparametros$.config)) %>% 
  map(~ Wokflows_set_map %>% 
        collect_predictions() %>% 
        filter(
          str_detect(wflow_id, .x[[1]]) & 
            .config== .x[[2]]
        )) %>% 
  
  map(~ .x %>%
        mutate(Model= str_extract(wflow_id, '^.{2}') ) %>% 
        select(Model,Var_respuesta, .pred_0, .pred_1, .pred_class)
  ) %>% 
  set_names(list('Escala','Cognicion','Bioquimica','Popurri') ) %>%
  map(~ .x %>% roc_curve(truth = Var_respuesta, .pred_0)) %>% 
  map2(., list('Escala','Cognicion','Bioquimica','Popurri'),
       ~ .x %>% mutate(model= .y)) %>% 
  bind_rows() 
  
Plot_ROC_CURVE_training_modelos_finales <- ROC_CURVE_training_modelos_finales %>% 
  ggplot(aes( x = 1 - specificity ,y = sensitivity, color = model ) ) +
  geom_path() +
  geom_abline(intercept = 0, slope = 1, linetype = 3) +
  theme_bw()


ROC_CURVE_testing_modelos_finales <- list(
  Modelo_Escalas_RandomForest_final,
  Modelo_Cognicion_RandomForest_final,
  Modelo_Bioquimica_RandomForest_final,
  Modelo_Popurri_RandomForest_final
  ) %>% 
  map(~.x %>% 
        collect_predictions() %>% 
        roc_curve(Var_respuesta, .pred_0)) %>% 
  map2(., list('Escala','Cognicion','Bioquimica','Popurri'),
       ~ .x %>%  mutate(model = .y)) %>% 
  bind_rows()

Plot_ROC_CURVE_testing_modelos_finales <- ROC_CURVE_testing_modelos_finales %>% 
  ggplot(aes(x = 1 - specificity ,y = sensitivity, color = model ) ) +
  geom_path() +
  geom_abline(intercept = 0, slope = 1, linetype = 3) +
  theme_bw()

ggpubr::ggarrange( Plot_ROC_CURVE_training_todos_los_posibles_modelos,
                   Plot_ROC_CURVE_training_modelos_finales, 
                   Plot_ROC_CURVE_testing_modelos_finales ,
                   nrow = 1 , 
                   common.legend = T,
                   legend = "bottom",
                   labels = c('Todos_training','Training','Test'))

14 Haciendo el machine learning entendible: SHAP Values

Los valores de Shapley (shap values), que deben su nombre al premio Nobel Lloyd Shapley, son un concepto de la teoría de juegos cooperativos que ha acabado encontrado aplicaciones en el aprendizaje automático.

En teoría de juegos, los valores de Shapley son la distribución equitativa de la contribución de cada jugador en un juego cooperativo entre todos los jugadores. En el contexto del aprendizaje automático, las características o variables se consideran jugadores, y el juego consiste en cuánto contribuye cada característica a la predicción de un modelo.

Es decir, los valores de Shapley son sencillamente las marginales de cada variable dentro del modelo. Compilarlos es computacionalmente extenso, ya que hay que hacer cada combinación de valores para una fila y evaluar de cuanto es su marginal. El valor final llamado Shap value es promedio de todas las aportaciones de una variable marginal.

Aparte de su extenso coste computacional, todo shap value asume que está utilizando un modelo de caja negra. de manera que si el modelo no está bien especificado en tema de interacciones, etc. es posible que estén dando indicaciones erróneas.

Actualmente se pueden compilar en tidymodels utilizando la librería agua, al ser una extensión de tidymodels del entorno de h2o. El problema es que esto hace que el motor de los modelos deba ser a la fuerza h2o.

Una propuesta sería buscar un modelo bien calibrado con tidymodels y pasarlo por h2o para hacer un diagnostico más eficaz.

h2o_start()
h2o::h2o.no_progress()

modelo_h20 <- rand_forest() %>% 
  set_engine("h2o", max_runtime_secs = 20) %>% 
  set_mode('classification') %>% 
  fit(Var_respuesta ~ ., data = Receta_modelo_Bioquimica %>% bake(new_data = NULL) )

entrada_h20 <- h2o::as.h2o(
  Receta_modelo_Bioquimica %>% bake(new_data = NULL) %>%
    dplyr::rename(.outcome = Var_respuesta))

modelo_h20 %>%  
  extract_fit_engine() %>% 
  h2o::h2o.shap_summary_plot(entrada_h20)

# modelo_h20 %>%  
#   extract_fit_engine() %>% 
#   h2o::h2o.explain(entrada_h20)

h2o_end()

Para cada predicción, los valores de Shapley representan la contribución de cada característica a la diferencia entre la predicción del modelo y la predicción media. En otras palabras, indican cuánto añade o resta cada característica a la predicción media.

Si una predicción media es de 0.5, y el shap value es de 0.1, indicará que esa variable sola podría llegar a tener un impacto tal que haga que la predicción llegue a 0.6.

Taller de tidymodels

Jornadas de R

Aitor Gonzalez, estadístico y data scientist en el dpt de medicina de la UAB

1 Antes de empezar, conozcámonos

2 ¿Por que tidymodels?

2.1 Un poquito de la historia

2.2 Un breve ejemplo

2.3 ¿Que aporta de tidymodels?

3 Presentando los datos

4 ¿De qué se compone tidymodels?

5 Preparando la partición de datos con Rsampling

5.1 Estratificar

5.2 Deshacer la partición

5.3 folds, creando conjuntos de validación

6 Pre-procesamiento de datos con `parsnip`

6.1 ¿Cómo cocinamos una receta?

6.2 Un par de ejemplos más

6.2.1 Imputación

6.2.2 PCA

6.2.3 Re-remuestrear

6.3 Caso práctico: receta

7 ¿Dónde está el modelo? Creando modelos con parnsip

7.1 Caso práctico, modelo de Random Forest

8 Unificando todo: Workflow

9 Eligiendo métricas con `yardstick`

10 Eligiendo hiperparámetros: tune y dials

10.1 Montando un espacio de híper-parametros: familia grid_

10.2 Ajustando los hiper-parametros: convergen `yardstick` con `dial` y `tune`

11 La conversión de todo, `fit` y `predict`

11.1 `last_fit`

12 Rendimiento del umbral (threshold performance en modelos de clasificación)

13 El cielo es el límite, `workflow_set()`

13.1 Todas las recetas

13.2 `workflow_set()` para unificarlas todas

13.3 Métricas de todos los modelos

13.4 Cruvas ROC todos los modelos

14 Haciendo el machine learning entendible: SHAP Values

Taller de tidymodels

Jornadas de R

Aitor Gonzalez, estadístico y data scientist en el dpt de medicina de la UAB

1 Antes de empezar, conozcámonos

2 ¿Por que tidymodels?

2.1 Un poquito de la historia

2.2 Un breve ejemplo

2.3 ¿Que aporta de tidymodels?

3 Presentando los datos

4 ¿De qué se compone tidymodels?

5 Preparando la partición de datos con Rsampling

5.1 Estratificar

5.2 Deshacer la partición

5.3 folds, creando conjuntos de validación

6 Pre-procesamiento de datos con parsnip

6.1 ¿Cómo cocinamos una receta?

6.2 Un par de ejemplos más

6.2.1 Imputación

6.2.2 PCA

6.2.3 Re-remuestrear

6.3 Caso práctico: receta

7 ¿Dónde está el modelo? Creando modelos con parnsip

7.1 Caso práctico, modelo de Random Forest

8 Unificando todo: Workflow

9 Eligiendo métricas con yardstick

10 Eligiendo hiperparámetros: tune y dials

10.1 Montando un espacio de híper-parametros: familia grid_

10.2 Ajustando los hiper-parametros: convergen yardstick con dial y tune

11 La conversión de todo, fit y predict

11.1 last_fit

12 Rendimiento del umbral (threshold performance en modelos de clasificación)

13 El cielo es el límite, workflow_set()

13.1 Todas las recetas

13.2 workflow_set() para unificarlas todas

13.3 Métricas de todos los modelos

13.4 Cruvas ROC todos los modelos

14 Haciendo el machine learning entendible: SHAP Values

6 Pre-procesamiento de datos con `parsnip`

9 Eligiendo métricas con `yardstick`

10.2 Ajustando los hiper-parametros: convergen `yardstick` con `dial` y `tune`

11 La conversión de todo, `fit` y `predict`

11.1 `last_fit`

13 El cielo es el límite, `workflow_set()`

13.2 `workflow_set()` para unificarlas todas